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ABSTRACT 

A large amount of observations have constrained cosmological parameters and the initial density 
fluctuation spectrum to a very high accuracy. However, cosmological parameters change with time 
and the power index of the power spectrum varies with mass scale dramatically in the so-called 
concordance ACDM cosmology. Thus, any successful model for its structural evolution should work 
well simultaneously for various cosmological models and different power spectra. 

We use a large set of high-resolution A^-body simulations of a variety of structure formation models 
(scale-free, standard CDM, open CDM, and ACDM) to study the mass accretion histories, the mass 
and redshift dependence of concentrations and the concentration evolution histories of dark matter 
halos. We find that there is significant disagreement between the much-used empirical models in the 
literature and our simulations. Based on our simulation results, we find that the mass accretion rate of 
a halo is tightly correlated with a simple function of its mass, the redshift, parameters of the cosmology 
and of the initial density fluctuation spectrum, which correctly disentangles the effects of all these 
factors and halo environments. We also find that the concentration of a halo is strongly correlated 
with the universe age when its progenitor on the mass accretion history first reaches 4% of its current 
mass. According to these correlations, we develop new empirical models for both the mass accretion 
histories and the concentration evolution histories of dark matter halos, and the latter can also be 
used to predict the mass and redshift dependence of halo concentrations. These models are accurate 
and universal: the same set of model parameters works well for different cosmological models and for 
halos of different masses at different redshifts, and in the ACDM case the model predictions match 
the simulation results very well even though halo mass is traced to about 0.0005 times the final mass, 
when cosmological parameters and the power index of the initial density fluctuation spectrum have 
changed dramatically. Our model predictions also match the PINOCCHIO mass accretion histories 
very well, which are much independent of our numerical simulations and of our definitions of halo 
merger trees. These models are also simple and easy to implement, making them very useful in 
modeling the growth and structure of dark matter halos. We provide appendices describing the step- 
by-step implementation of our models. A calculator which allows one to interactively generate data for 
any given cosmological model is provided at http : //www . shao . ac . cn/ dhzhao/mandc . html, together 
with a user-friendly code to make the relevant calculations and some tables listing the expected 
concentration as a function of halo mass and redshift in several popular cosmological models. We 
explain why ACDM and open CDM halos on nearly all mass scales show two distinct phases in their 
mass growth histories. We discuss implications of the universal relations we find in connection to the 
formation of dark matter halos in the cosmic density field. 

Subject headings: cosmology: miscellaneous — galaxies: clusters: general — methods: numerical 



1. INTRODUCTION 

In the current cold dark matter (CDM) paradigm of 
structure formation, a key concept in the buildup of 
structure in the universe is the formation of dark mat- 
ter halos. These halos are quasi-equilibrium systems of 
CDM particles formed through nonlinear gravitational 
collapse in the cosmic density field. Since galaxies and 
other luminous objects are believed to form by cooling 
and condensation of the baryonic gas in potential wells 
of dark matter halos, the understanding of the forma- 
tion and properties of CDM halos is an important part 
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of galaxy formation. 

One of the most important properties of the halo pop- 
ulation is their density profiles. Based on A'^-body simu- 
lations, Navarro, Frenk, & White (1997, hereafter NEW) 
found that CDM halos can in general be approximated 
by a two-parameter profile. 
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where Rs is a characteristic "inner" radius at which the 
logarithmic density slope is —2, and ps is the density at 
Rg. A halo is often defined so that the mean density 
Ph within the halo radius Rh is a factor A/i times the 
mean density of the universe (p) at the redshift (z) in 
consideration. The halo mass can then be written as 
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One can define the circular velocity of a halo as 
{GMh/RhY''^, so that 
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where H{z) is the Hubble parameter and ri„j(z) is the 
cosmic mass density parameter at rcdshift z. The shape 
of an NFW profile is usually characterized by a concen- 
tration parameter c, defined as c = Rh/Rs- It is then 
easy to show that 
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We denote the mass within Rg by Mg, and the circular 
velocity at Rg by Vg- These quantities are related to c 
and Mh by 
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In the literature, a number of definitions have been 
used for Ah (hence Rh)- Some authors opt to use 
= 200 (e.g., Jenkins et al. 2001) or = 200/f]„(z) 
(e.g., NFW), while others (e.g.. Bullock et al. 2001; Jing 
& Suto 2002) choose Ah = Avir according to the spheri- 
cal virialization criterion (Kitayama & Suto 1996; Bryan 
& Norman 1998). These different definitions can lead to 
sizable differences in c for a given halo, and the differ- 
ences are cosmology-dependent. In our discussion in the 
main text we use Ah = Avir, and adopt the fitting for- 
mulae of Avir obtained by Bryan & Norman (1998). The 
value of Avir ranges from ~ 180 at high redshift to ~ 340 
at present for the current 'concordance' ACDM cosmol- 
ogy. In appendices A and B, we also provide results 
based on the other two definitions of Ah- To avoid con- 
fusing, we denote halo mass, radius, concentration and 
circular velocity by M, i?, c, and V for the last definition 
{Ah = Avir), by M200C, -R200C, C200C, and ¥2000 for the sec- 
ond definition [Ah = 2QQ/Q.^{z)), and by M2ooin, -R200m, 
C200m, and V200m for the first definition {Ah = 200). 

The structure of a halo is expected to depend not only 
on cosmology and fluctuation power spectrum, but also 
on its formation history. Attempts have therefore been 
made to relate the halo concentration to quantities that 
characterize the formation of the halo. In their original 
paper, NFW suggested that the characteristic density 
of a halo, pg, is a constant (fc) times the mean den- 
sity of the universe, p{zf), at the redshift Zf (referred 
to as the formation time of the halo by NFW) where 
half of the halo's mass was first assembled into progen- 
itors more massive than / times the halo mass. NFW 
used the extended Press-Schechter formula to calculate 
Zf and found that the anticorrelation between c and M 
observed in their simulation can be reproduced reason- 
ably well with a proper choice of values of the parameters 
k and /. 

Subsequent investigations demonstrated that addi- 
tional complexities are involved in the CDM halo struc- 
ture. First, halos of a fixed mass may have significant 
scatter in their c values (e.g., Jing 2000), although there 
is a mean trend of c with M. If this trend is indeed 
due to a correlation between concentration and forma- 
tion time, the scatter in c may reflect the expected scat- 



ter in the formation history for halos of a given mass 
(e.g., Jing 2000, Lu et al. ' 2006). Second, Bullock et 
al. (2001, hereafter BOl) and Eke et al. (2001, hereafter 
EOl) found that the halo concentration at a flxed mass 
is systematically lower at higher redshift. BOl proposed 
a model with c oc {M/M^,)""-'^^ /{I + z) (M* being the 
mass at which the rms of the linear density field is ^ 1), 
which has a stronger halo-mass dependence, and a much 
stronger redshift dependence, than that predicted by the 
NFW model. EOl proposed a similar model that has a 
weaker mass dependence, and a slightly weaker redshift 
dependence, than the BOl model. Both of these models 
have been widely used in the literature to predict the 
concentration halo mass relation. Using the same set of 
numerical simulations as that used in BOl, Wechsler et 
al. (2002; hereafter W02) found that, over a large mass 
range, the mass accretion histories (hereafter MAHs) of 
individual halos identified at redshift ^obs can be approx- 
imated by a one-parameter exponential form, 

M{z) = M(2obs)exp[-2(z - Zobs)/(l + ^/)] , (6) 

where -2/ is the formation time of a halo, determined 
by fitting the simulated MAH with the above formula. 
Assuming that c equals 4.1 at Zf and grows proportion- 
ally to the scale factor a cx 1/(1 + z), W02 proposed 
a recipe to predict c for individual halos through their 
MAHs, and found that their model can reproduce the 
dependence of c on both halo mass and redshift found 
in BOl. Using a large set of Monte Carlo realizations of 
the extended Press-Schechter formalism (EPS), van den 
Bosch (2002) showed that the average MAH of CDM 
halos follows a simple universal function with two scal- 
ing variables. Based on a set of high-resolution A^-body 
simulations, Zhao et al. (2003a, 2003b, hereafter Z03a 
and Z03b, respectively) found that the MAH of a halo 
in general consists of two distinct phases: an early fast 
phase and a late slow phase (see also Li et al. 2007 
and Hoffman et al. 2007). As shown in Z03a, the early 
phase is characterized by rapid halo growth dominated by 
major mergers, which effectively reconfigure and deepen 
the gravitational potential wells and cause the collision- 
less dark matter particles to undergo dynamical relax- 
ation and mix up sufficiently to form the core structure, 
while the late phase is characterized by slower quiescent 
growth predominantly through accretion of material onto 
the halo outskirt, little affecting the inner structure and 
potential. Z03a proposed that the concentration evolu- 
tion of a halo depends much on its mass accretion rate 
and the faster the mass grows, the slower the concen- 
tration increases. In particular, they predicted that ha- 
los which are still in the fast growth regime should all 
have a similar concentration, c ~ 4. Using a combina- 
tion of A^-body simulations of different resolutions, Z03b 
studied in detail how the concentrations of CDM halos 
depend on halo mass at different redshifts. They con- 
firmed that halo concentration at the present time de- 
pends strongly on halo mass, but their results also show 
marked differences from the predictions of the BOl and 
EOl models. The mass dependence of halo concentra- 
tions becomes weaker at higher redshifts, and at z > 3 
halos with M > 10^^h~^M.Q all have a similar median 
concentration, c ~ 4. While the median concentrations 
of low-mass halos grow significantly with time, those of 
massive halos change very little with redshift. These re- 
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suits arc in good agreement with the empirical model 
proposed by Z03a and favored by the Chandra observa- 
tion of Schmidt & Allen (2007), but are very different 
from the predictions of the models proposed by BOl, EOl 
and W02. Recently, Gao et al. (2008) and Maccio et al. 
(2008) confirmed the results of Z03b with the use of a 
large set of simulated halos, suggesting again that the 
much-used BOl and EOl models need to be revised. 

There have been attempts to modify the early em- 
pirical models of halo MAH and halo concentration to 
accommodate the new simulation results. Miller et al. 
(2006) and Neistcin et al. (2006) provided MAH mod- 
els based on theoretical EPS formalism. With help of 
the Millennium simulation, Ncistein & Dekel (2008) pre- 
sented a empirical algorithm for constructing halo merger 
trees, which is significantly better than their previous 
method based on EPS but fails to reproduce the non- 
Markov features of merger trees. Moreover, because 
of the limited resolution of the simulation, this algo- 
rithm is tmicd to reproduce MAHs within redshift 2.5 
for only massive halos. McBride et al. (2009) examined 
halo MAHs of the Millennium sirmilation and modeled 
the scatter among different halos by fitting individual 
MAHs with a two-parameter function, which was proved 
to be more accurate than the one-parameter function 
of W02. Note that both the two models above are re- 
stricted to the specific cosmology and the specific fluc- 
tuation power spectrum of the Millennium simulation, 
similar to WMAPl results. Gao et al. (2008) confirmed 
the finding of BOl that the NFW prescription overpre- 
dicts the halo concentrations at high redshift, and tried 
to overcome this shortcoming by modifying the defini- 
tion of halo formation time. However, the revised model 
still fails to reproduce the evolution of concentrations of 
galaxy sized halos. Modifying the parameters of the BOl 
and EOl models can also reduce the discrepancy between 
these models with simulation results at high masses, but 
the overly rapid redshift evolution of the concentration 
remains. Motivated partly by Z03a and Z03b, Maccio 
et al. (2008) presented a model based on some modifi- 
cation of the BOl model, and found that it can repro- 
duce the concentration-mass relation at z = in their 
simulations. Unfortunately, their model is not universal, 
because the normalization of the concentration-mass re- 
lation has to be calibrated for each cosmological model 
and for each redshift and thus the model cannot be used 
to predict correctly the redshift evolution of halo concen- 
trations. In addition, as we will show below, the model 
of Maccio et al. has the same shortcoming as the original 
BOl model in that it predicts too steep a concentration- 
mass relation for halos at high redshift in the current 
ACDM models and for halos at all redshift in the 'stan- 
dard' CDM model with il„i.o = 1- Thus, all the exist- 
ing empirical models for the halo concentration can at 
best provide reliable predictions for dark matter halos in 
limited ranges of halo mass and redshift (often around 
^; = 0), and only for some specific cosmological models 
(according to which the empirical models are calibrated). 
Clearly, they are insufficient in the era of precision cos- 
mology. 

We believe that even though a large amount of observa- 
tions have constrained cosmological parameters and the 
initial density fiuctuation spectrum to a very high ac- 
curacy, cosmological parameters change with time and 



the power index of the power spectrum varies with mass 
scale dramatically in the so-called concordance ACDM 
cosmology and thus any successful model for its struc- 
tural evolution should work well simultaneously for var- 
ious cosmological models and different power spectra. 

In this paper, we use A/'-body simulations of a variety 
of structure formation models, including a set of scale- 
free (SF) models, two ACDM (LCDM) models, the stan- 
dard CDM (SCDM) model, and an open CDM (OCDM) 
model, to investigate in detail the MAHs, the mass and 
redshift dependence of concentrations and the concentra- 
tion evolution histories of dark matter halos. We show 
that early empirical models for halo MAHs and concen- 
trations all fail significantly to describe the simulation 
results. Based on our simulation results, we develop new 
empirical models for both the MAHs and the concen- 
tration evolution histories of dark matter halos, and the 
latter can also be used to predict the mass and redshift 
dependence of halo concentrations. These models are ac- 
curate and universal, in the sense that the same set of 
model parameters works well for different cosmological 
models and for halos of different masses and at differ- 
ent redshifts. These models are also simple and easy to 
implement, making them very useful in modeling the for- 
mation and structure of dark matter halos. Furthermore, 
the universal relations we find may also provide impor- 
tant insight into the formation processes of dark matter 
halos in the cosmic density field. 

The organization of the paper is as follows. In 
Section 2, we present the set of A/'-body simulations 
used in the paper. We describe the simulated MAHs 
and the corresponding modeling in Section 3. Our 
modeling of halo concentrations is presented in Sec- 
tion 4. Finally, we discuss and summarize our re- 
sults in Section 5. We provide appendices describ- 
ing the step- by-step implementation of our models. 
A calculator which allows one to interactively gen- 
erate data for any given cosmological model is pro- 
vided at http : //www. shao .ac . cn/dhzhao/mcindc .html, 
together with a user-friendly code to make the relevant 
calculations and some tables listing the expected concen- 
tration as a function of halo mass and redshift in several 
currently popular models of structure formation. 

2. NUMERICAL SIMULATIONS AND DARK MATTER 
HALOS 

We use a very large set of cosmological simulations 

to study the formation and structure of dark matter 
halos of different masses in different cosmological mod- 
els. The first subset contains the simulations that were 
used in Z03b: the then 'concordance' ACDM model with 
density parameter f2m,o = 0.3 and a cosmological con- 
stant given by Oa,o = 0.7. These simulations, labeled 
as LCI)Ml-3 in Table 1, were generated with a parallel- 
vectorized Particle-Particle/Particle-Mesh code (Jing & 
Suto 2002). The linear power spectrum has a shape pa- 
rameter F = flm,oh = 0.20, and an amplitude specified 
by (Tg = 0.9, where h is the Hubble constant in units of 
100kms~^Mpc^^, and is the rms of the linear den- 
sity field smoothed within a sphere of radius 8ft-~^Mpc 
at the present time. We used 256'^ particles for the sim- 
ulation of boxsize L — 25/i^^Mpc, and 512'^ particles 
for the other two simulations, L = 100 and 300/i~^Mpc 
(see Table 1). The simulations with L = 25 h~^Mpc and 
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Fig. 1. — rms of the linear density fluctuation a-{AI) as a function 
of the mass scale M for scale-free simulations with n = 1, 0, — 1, 
and -2, and for ACDM simulations LCDMl-3. 

100 h^^Mpc were evolved by 5000 time steps with a force 
softening length rj (the diameter of the S2 shaped parti- 
cles, Hockney & Eastwood 1981) equal to 2.5ft.^^kpc and 
10/i~^kpc, respectively. In order to examine the depen- 
dence on cosmological parameters, we also use simula- 
tions for the 'Standard' CDM model (rim.o — 1) and the 
open CDM model (r^m.o — 0.3, r^A.o = 0), as described 
in Jing & Suto (2002, 1998). These are listed as SCDMl, 
SCDM2 and OCDMl in Table 1, together with the corre- 
sponding model and simulation parameters. In addition, 
we also use a set of four scale-free (SF) simulations in 
an Einstein de Sitter universe, with the linear power-law 
power spectra given by P{k) (x k". These models are 
listed as SFl-4 in Table 1 and Figure 1 presents the rms 
of the linear density field a{M) as a function of the mass 
scale M for SFl-4 and LCDMl-3. 

We use both the friends of friends (hereafter FOE) al- 
gorithm and the spherical overdensity (hereafter SO) al- 
gorithm of Lacey & Cole (1994) to identify dark matter 
groups in the simulations. The EOF algorithm connects 
every pair of particles closer than 0.2 times the mean 
particle distance and the SO algorithm selects spherical 
regions whose average density is equal to Avir times the 
mean cosmic density p. We select groups at a total of 20- 
30 snapshots for each CDM simulation and at 10 snap- 
shots for each SF simulation. The outputs are uniformly 
spaced in the logarithm of the cosmic scale factor a. For 
each group, we choose the particle of the highest local 
density as the group center, and around it we select as 
'halo' a spherical region whose average density is AvirP, 
with a routine very similar to the SO algorithm except 
that here the group center is fixed. We find that the re- 
sults presented in this paper do not depend significantly 
on the group finding algorithm. Our following presen- 
tation is based on the EOF groups and we use all EOF 
groups without applying any further selection criteria. 

The simulations described above are also used to study 
the density profiles and the concentrations of dark mat- 
ter halos. We use all halos containing more than 500 
particles for our analysis. Each halo is fitted with the 
NFW profile, using a similar method as described in 
Z03b. As demonstrated in Z03b, 500 particles are suffi- 



10 



0.1 

0.01 
10'" 

1015 



^ I 


1 




: 


\ 

x\ SFl 




r 


i-W SF2 


-_ 


- 




- 








- 




- 


; 












H 


tj \ — 1 1 1 1 II 

LCDMl- 


3 - 




SCDMl- 


21 


r 




















1 





z+1 

Fig. 2. — Median halo MAHs obtained from scale-free simulations 
with n = 1, 0, -1, and -2 and from CDM simulations LCDMl-3, 
SCDMl— 2, and OCDMl. Different lines represent median histories 
with different final halo masses. For each final mass, the result is 
plotted to a redshift when the histories of 90% of the halos can still 
be traced. 

cient for measuring the concentration parameter reliably. 
In order to examine concentrations of very massive halos, 
we also use three additional simulations with large box- 
sizes (see Jing, Suto & Mo 2007). These simulations are 
listed as LCDM4-6 in Table 1. Note that the model pa- 
rameters of these simulations are slightly different from 
those of LCDMl-3. Furthermore, the initial power spec- 
trum for these simulations is obtained using the fitting 
formula of Eisenstein & Hu (1998), which includes the 
effects of baryonic oscillations. The details of these sim- 
ulations can be found in Jing, Suto & Mo (2007), and 
all the change in model parameters are properly taken 
into account when we compare our model predictions of 
halo concentrations with the simulation results. For the 
same reason, we also use a SCDM simulation of boxsize 
L = 300/i-iMpc (listed as SCDM3 in Table 1) to in- 
vestigate the concentration-mass relation at the highest 
halo mass end in this model. 

3. A UNIVERSAL MODEL FOR THE HALO MASS 
ACCRETION HISTORIES 

Since we want to study how a halo grows with time, 
we need to construct the main branch of the merger tree 
for each halo. Given a group of dark matter particles at 
a given output time (which we refer to as group 2), we 
trace all its particles back to an earlier output time. A 
group (group 1) at the earlier output is selected as the 
"main progenitor" of group 2 if it contributes the largest 
number of particles to group 2 among all groups at this 
earlier output. We found that in most cases more than 
half of the particles of group 1 is contained in group 2. 
We refer to group 2 as the "main offspring" of group 1. 
We construct the main branch of the merger tree, for each 
of the most massive 10000 halos identified at redshift 
z — 0, until the number of particles in the progenitor 
drops below 10 or the main progenitor cannot be found 
anymore. More than 97% of the histories analyzed here 
can be traced until the particle number goes below 10. 

As an illustration. Figure 2 shows the median MAHs of 
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TABLE 1 

A Summary of Simulation Parameters 



Simulation 




h 


0"8 




f^A,0 


Np 


L[h-'-Mpc] 


Nnq 


ri 




SFl 


1 




1 


1 





512-^ 


160 


512 


0.016 


1481.3 


SF2 







1 


1 





512^ 


160 


512 


0.016 


363.2 


SF3 


-1 




1 


1 





512=* 


160 


512 


0.016 


75.0 


SF4 


-2 




1 


1 





512=* 


160 


512 


0.016 


11.3 


LCDMl 


1 


0.667 


0.9 


0.3 


0.7 


256=* 


25 


512 


0.0025 


72 


LCDM2 


1 


0.667 


0.9 


0.3 


0.7 


5123 


100 


512 


0.01 


72 


LCDM3 


1 


0.667 


0.9 


0.3 


0.7 


512=* 


300 


512 


0.03 


36 


LCDM4 


1 


0.71 


0.85 


0.268 


0.732 


1024^ 


300 


1024 


0.01 


72 


LCDM5 


1 


0.71 


0.85 


0.268 


0.732 


10243 


1200 


1024 


0.072 


72 


LCDM6 


1 


0.71 


0.85 


0.268 


0.732 


10243 


1800 


1024 


0.144 


72 


OCDMl 


1 


0.667 


1.0 


0.3 





2563 


50 


256 


0.02 


72 


SCDMl 


1 


0.5 


0.55 


1 





2563 


25 


512 


0.0025 


72 


SCDM2 


1 


0.5 


0.55 


1 





5123 


100 


256 


0.02 


72 


SCDM3 


1 


0.5 


0.55 


1 





5123 


300 


256 


0.03 


36 



Note. — The Bardeen et al. (1986) fitting formula for transfer function is used for all CDM simulations except LCDM4-6. For LCDM4- 
6, the fitting formula of Eisenstein & Hu (1998) is adopted, with baryon density parameter SI5 q = 0.045. Note that a short-wavelength 
cutoff in the linear power spectrum must be included in the SFl simulation. 



halos both in the SF simulations and in the CDM models. 
The halo mass is normalized by the characteristic mass 
Af^ the final output in the SF simulations. Clearly, 
MAHs depend on the power spectrum. This is because 
halos grow faster in a model with a smaller value of n. 
The effect of cosmological parameters can also be seen 
clearly in CDM models: the growth of halos is slower in 
the LCDM model than in the SCDM model (despite the 
fact that the LCDM has more clustering power than the 
SCDM) and is the slowest in the OCDM. In our analysis, 
we have taken into account the incompleteness effect of 
the MAHs at the earliest epochs when the progenitors 
containing less than 10 particles cannot be followed in 
the simulation. This effect leads to an overestimation of 
MAHs. In order to correct for it, statistical analysis of 
the MAHs is carried out only up to the redshift where 
the progenitors of 90% of all the halos in consideration 
can be traced. 

As one can see from Figure 2, the MAHs show a com- 
plex dependence on power spectrum, cosmology and halo 
mass. The goal of this section is to find an order in such 
complexity so as to obtain a universal model to describe 
aU the MAHs. 

3.1. The Expected Asymptotic Behavior at High Redshift 

Before modeling the MAHs in detail, let us first con- 
sider some generic properties of the growth of CDM halos 
in the cosmic density field. In the CDM scenario, struc- 
tures form in a hierarchical fashion, and the growth of 
dark matter halos in general has the following properties. 
First, at any given time, more massive halos are, on av- 
erage, growing faster because they sit in higher density 
environments. Second, the higher the redshift, the faster 
the halos of a fixed mass grow because at high redshift 
they are relatively more massive with respect to others. 
Third, the MAH of a halo depends on power spectrum 
of the initial density fluctuation, as mentioned before. 
Fourth, the MAH of a halo also depends on cosmology, 
because of the change of the background expansion. The 
last, halo growth also suffers from some nonlinear pro- 
cesses in their local environments. All these factors en- 
tangle together and keep varying in their own compli- 
cated ways. For example, in a ACDM universe, the cos- 
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Fig. 3. — Evolution of cosmological density parameter Qrn and 
cosmological constant parameter Q\ for a flat universe with cos- 
mological constant. 



mological density parameter f2„i decrease with time and 
the power index of the initial density fluctuation spec- 
trum increase with mass scale dramatically, as shown in 
Figure 3 and Figure 1, respectively. This is why one 
need a universal evolution model which works well for 
various cosmologies and different spectra, as mentioned 
in the introduction already. For exactly the same rea- 
son, however, we may not expect that variations of all 
these factors compensate exactly to bring us a universal 
function form for MAH in terms of M and z. Thus, we 
turn to model halo mass accretion rate instead, which 
can be integrated to build the MAH, and in order to de- 
scribe the accretion rate accurately and universally, we 
need to disentangle all effects due to halo mass, redshift, 
the cosmological parameters, the linear power spectrum 
and those relevant nonlinear processes. 

Consider two average halos of different masses at a 
given time. Since the more massive one grows faster, 
we expect that at a slightly earlier time, the two halos 
were closer in mass (in terms of the mass ratio), and 
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the masses of their progenitors at higher redshift are ex- 
pected to be even closer. The two halos are expected to 
have similar mass accretion rate at sufBciently high red- 
shift, because otherwise with finite difference in accretion 
rates their mass accretion trajectories would cross. This 
would lead to an unlikely situation where more massive 
progenitors actually grow slower than low-mass progen- 
itors at redshifts higher than the crossing redshift. This 
means that the average M AHs of halos of different masses 
should all have the same asymptotic behavior at high 
redshift. Here we try to look for this behavior. 

Let us first consider an Einstein de Sitter universe with 
an initial Gaussian density fluctuation field of a power- 
law power spectrum, P{k) cx k". As there is no char- 
acteristic scale either in time or in space, this is a case 
where the structure formation is a self-similar process. 
When halo mass is scaled with a time-dependent charac- 
teristic mass, all statistical quantities as functions of the 
scaled mass should be the same at all redshifts. Given 
a linear power spectrum, we can linearly extrapolate to 
z = 0, and estimate, the rms of the linear density field 
a{M) on a mass scale M. According to the spherical col- 
lapse model, the linear critical overdensity for collapse 
at redshift z is 6c{z) = 6'^[Qm{z),V,A{z)]/ D{z), where 
D{z) is the linear growth factor, which is 1/(1 -I- z) in an 
Einstein de Sitter universe, and 5'^[^lm{z),i^A{z)] is the 
conventional critical overdensity for collapse at redshift 
z, which is a constant ^ 1.686 in this universe. Because 
Sc{z) is a function of only z and cr(M) is a function of only 
M, we can think of l/a'^{M) as the "mass" and l/Sl{z) 
as the "time". Furthermore, since M^, which satisfies 
a{Mi,) = Sc{z), represents the characteristic nonlinear 
mass at redshift z, 1/S'^{z) actually corresponds to the 
characteristic nonlinear "mass" at redshift z if 1/(T^(M) 
is regarded as "mass" . In this case, the scaled "mass" is 

= 5l{z) / a"^ {M) , i.e., the peak height of the halo in the 
conventional terminology. 

Suppose that there are some halos at a given 
redshift Zobs and the average mass growth rate is 
d\na^{M)/d\n5l{z%^^^ = 1, i.e., d\nv / d\n5l{z)l^^^ = 
so that V — constant. Then at a slightly higher redshift, 
zi , the scaled mass v of their progenitors will be the same 
as that at Zobs- Since the statistical properties of halos 
of the same scaled mass are the same for all redshifts in 
the self-similar case, these progenitors will have the same 
growth rate as that at Zobs, i-e-, both the mass growth 
rate d\na^{M)/d\nS-^{z)\z^^^^, and the scaled mass v, re- 
mains constant with time. Therefore, the main progeni- 
tors at different redshifts will all have the same v and the 
same growth rate as they have at redshift Zobs, and the 
average MAH will be a straight line in the logarithmic di- 
agram of a versus Sc- This line is exactly the asymptotic 
behavior we are looking for. 

In the case of cold dark matter cosmology (SCDM, 
LCDM, and OCDM), structure formation is not self- 
similar. Nevertheless, we can still use l/a^(M) to rep- 
resent halo mass, l/5'^{z) to represent time, and v to 
represent the scaled mass. However, it is not guaranteed 
that the a-Sc relation has a unique asymptotic behavior 
as it has in the SF case. As we will see below, simu- 
lated median MAHs of different final masses show the 
same asymptotic behavior at high z not only for the SF 
models but also for the CDM models. We can then use 
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Fig. 4. — Median (solid blue) and scatter (dashed green) of the 
simulated MAHs, in terms of l/a{M) vs. Sc{z), for halos of three 
final masses in the scale-free case with n = — 1. Dashed green lines 
enclose 60% of the histories in each mass bin. 

the result to guide our modeling of the halo MAHs in 
simulations. 

3.2. Toward a Universal Model for the Halo Mass 
Accretion Histories 

If we use l/a{M) as mass and l/6c{z) as time, the 
MAH looks like the curves shown in Figure 4, where 
we have adopted the n — —1 SF model as an exam- 
ple. The figure clearly shows the asymptotic behavior 
described above. The MAHs of halos of different present 
masses converge at high redshift, and the scaled mass 
= Sc{z)/(7{M) approaches a constant, as one can see 
in the case of the most massive halos in the figure, and as 
expected from the asymptotic behavior discussed above. 

Figure 5 shows the median accretion rate, 
da{M)/dSc{z), as a function of the scaled mass, 
y/u = Sc{z) /a-{M), for the four SF models. For each 
model, a set of symbols of the same type connected 
with a solid line represents accretion rates of progenitors 
of the halos of a given mass at z = 0, and the filled 
symbol represents the end of the MAH at z = 0. It is 
interesting to note that, for each SF model, the results 
for different halo masses are similar except at the end 
of the history. The accretion rate plotted in this way 
is approximately a constant for all the progenitors. As 
mentioned above, in each model the scaled mass is 
asymptotically a constant at high z. This corresponds 
to a line, da{M)/dSc{z) = a{M)/Sc{z), which is plotted 
as the dashed line in the figure together with a star 
showing the corresponding asymptotic accretion rate. 
At the end of each history, the mass accretion rate is 
somewhat lower than the average of their progenitors. 
This deviation from self-similarity is not a numerical ar- 
tifact, as the simulation resolution is sufficient at z = 0. 
Instead, it reflects the fact that the main progenitors are 
a special subset of the total halo population. While the 
halos of a given at the end of MAHs are chosen only 
by mass, their progenitors have an additional selection 
bias: they are chosen to be main progenitors of more 
massive halos. This biases the progenitor halos to have 
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Fig. 5. — Median mass accretion rate for our four scale-free mod- 
els. Purple, blue, dark blue and green symbols and lines represent 
results for scale-free simulations with n = 1, 0, — 1, and —2, respec- 
tively. Different symbols represent median histories with different 
final halo masses, and in each history we use a solid symbol to 
denote the final snapshot of the history. The four stars and the 
dashed curve represent the expected asymptotic behavior for high 
peaks in scale-free models. 

higher mass accretion rates than typical halos of a 
given V value. In other words, a fraction of halos in 
the total population are not main progenitors of larger 
halos in the future but instead will be swallowed by 
more massive halos. Their mass accretion rates may 
have been suppressed by their massive neighbors due to 
environment-heating or tidal-stripping, as envisaged in 
Wang et al. (2007), and thus their merger trees should 
also have some non-Markov features. 

As seen in Figure 5, accretion rates are different in dif- 
ferent SF models, again indicating that the accretion rate 
depends on the shape of the power spectrum. In order 
to model such dependence, we use d\!ia{M) / dh\5c{z) 
instead of da{M) / d6c{z) to represent the accretion rate. 
The results are shown in the upper left panel of Figure 
6, which again indicates that the accretion rate depends 
strongly on the power index of the linear spectrum. How- 
ever, we find that this dependence can be scaled away if 
we use the following variable instead of 5{z) / a{M): 

w{z,M) = S,{z)/s{M), (7) 

S{M) = a{M) X lOdlS<^/dlgrn\M (g) 

This is shown in the upper right panel of Figure 6. With 
the use of w{z, M), all the halo MAHs lie on top of each 
other, except for the snapshots close to the end of the 
MAHs. Even more remarkably, the accretion rates of 
the progenitors at high redshift can all be well described 
by a straight line. 



where 



d\gaiM)/d\gSc{z) 



w{z, M) 
5.85 



(9) 



which is shown in the panel as the solid line. Note that 
Ig s(M) is equal to \ga[M) plus its logarithmic slope at 
mass M: dlga / d\g7n\M- For a given power spectrum 
with n > —3, a{M) decreases with halo mass and so 
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Fig. 6. — Same as Figure 5, except that both the horizontal and 
vertical axes have been rescaled as indicated by the axis labels. 
Results for the six LCDM histories shown in Figure 2 with final 
halo mass larger than 10^'^ Mq are also plotted here, in the two 
lower panels. The black solid straight lines in the right two panels 
have a slope 1/5.85, and the dashed curve in each of the two panels 
represents the envelope given by halos in the end of the MAHs in 
simulations LCDMl— 3. 
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Fig. 7. — Same as Figure 6, except that the horizontal axis is 
rescaled again, as indicated by the axis label. Here p is defined in 
the text. 



w{z, M) increases with halo mass. The simple relation 
given by Equation (9) also describes well the MAHs at 
high redshifts in the LCDM model (see the lower two 
panels of Figure 6), although the power spectrum and 
the cosmology are very different from the SF models. 

As shown in the two right panels of Figure 6, the ac- 
cretion rate is systematically lower than that given by 
Equation (9) at the end of the MAHs. We find that this 
decrease of the accretion rate can be accounted for by 
replacing w{z, M) with w{z, M) — p{z, Zohs, Mobs)- The 
shift, p(z, Zobsj -Mobs), in the horizontal axis depends on 
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Fig. 8. — Predicted and simulated median halo MAHs for the 
scale-free model with n = — 1. Black symbols are simulation re- 
sults, and solid red lines represent our model predictions. Results 
are shown for halos in five bins of final mass. 

redshift and halo mass in the following way, 

p{z, Zobs, Mohs) = p{Zohs, Zohs, Mobs) 



xMax 



0,1- 



IgScjz) - IgScjZohs) 

0.272/w(zobs,Mobs) 



(10) 



where 



p(Zobs, ^obs, Mobs) = 



1 



U>(2:obs,-Mobs) 



l + [w(Zobs,Mobs)/4]f 



(11) 

is the shift at Zobs, the redshift at which the final halo is 
identified. The horizontal gap between the dashed and 
solid curves in the right two panels of Figure 6 shows 
this shift. As one can see, the shift describes well the 
deviation of the mass accretion rates in the end of the 
histories (the solid symbols) for the LCDM model. For 
the SF models, the shift appears to be too much, but 
this is mainly due to the fact that the interval of the last 
two snapshots in the SF models is quite big, and the ac- 
cretion rate is not estimated accurately. Figure 7 shows 
d\ga{M)/d\g6c{z) versus w{z, M) - p{z, Zohs, Mohs)- It 
is clear that the relation is much tighter, and is well de- 
scribed by the following relation. 



d\ga{M)/d\gSc{z) = 



W{z, M) - p{z, Zobs, Mobs) 

5.85 



(12) 



which is shown as the solid lines. The same results are 
also found for the SCDM and OCDM simulations, al- 
though they are not shown here. 

Given a cosmology and a linear power spectrum, it is 
easy to calculate s(M) for a halo of mass Mobs at a given 
redshift Zohs- From Equation (12), one can estimate the 
mass accretion rate, d\ga{M)/d\gSc{z) at Zobs- One can 
then compute cr(M) [or equivalently s{M)] at a redshift 
z incrementally higher than Zobs, thus tracing the MAH 
to higher redshifts. The MAH, in terms of s(M) and 
Sc{z), thus computed are shown in Figures 8 and 9. Since 
s(M) is a monotonic function of M and Sc{z) a mono- 
tonic function of z, the s{M)-6c{z) relation can easily be 




5o(z) 

Fig. 9. — Predicted and simulated median halo MAHs for simu- 
lations LCDMl— 3. Symbols and lines are the same as in Figure 8. 
Here results are shown for halos in seven bins of final mass. 

converted into an MAH in terms of AI and {1 + z). Fig- 
ures 10-12 show the median MAHs so obtained for halos 
of different final masses at different redshifts in various 
models (solid curves), in comparison with the simula- 
tion results (circles). For the ACDM case, the model 
predictions match the simulation results very well even 
though halo mass is traced to about 0.0005 times the final 
mass, when cosmological parameters and the power index 
of the initial density fiuctuation spectrum have changed 
dramatically. Moreover, the same set model parameters 
works pretty good for all halo masses, all redshifts and 
all cosmologies we are studying here. The typical errors 
of the predicted median MAHs in most cases are < 10%. 
Somewhat larger errors seen in the highest mass bin in 
Figure 11 may be due to the inaccurate determination 
in the simulation because the number of halos is small 
in this mass bin. These predictions are much more accu- 
rate than those of the model of W02, shown as the short- 
dashed lines in Figures 10-12,** and those of the model 
of van den Bosch (2002), shown as the long-dashed lines 
in Figures 10-12.^ Furthermore, unlike the latter two 
models, the MAHs predicted by our model for halos of 
different final masses do not cross at high redshift. Note 
that the W02 model works pretty well for low-mass halos 
in the ACDM and the OCDM models (Af < IO^^Mq). 
However, for more massive halos in these two models and 
for all halos in the scale-free and the SCDM simulations, 
it fails to provide an accurate description for the MAHs. 
Even though the W02 model and ours give similar MAHs 

* Following the instruction on Bullock's Web site 
http://www.physics.uci.edu/~bullock/CVIR/, we set F = 0.015 
in the BOl model (see Section 4.1 for more details) and use the 
returned collapsing redshift as input for the free parameter 2/ in 
the W02 model, Equation (6). 

^ Note that predictions of the model of van den Bosch (2002) 
are the average MAHs while the predictions of ours are the median 
MAHs. If the same definition is adopted, the difference between 
predictions of these two models will be even larger in the CDM 
cases for which our model always predicts more massive progeni- 
tors. This is because individual MAHs show a log-normal distri- 
bution (D. H. Zhao et al. 2010, in preparation) and so the average 
value should be somewhat higher than the median one. 
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Fig. 10. — Predicted and simulated median halo MAHs for scale- 
free models with n = 1, 0, —1 and —2. Black symbols are simula- 
tion results, while solid red lines represent our model predictions. 
Those histories ended at high redshifts are MAHs with final ha- 
los chosen at those redshifts. The short-dashed blue lines and the 
long-dashed green lines are predictions of the model of W02 and 
predictions of the model of van den Bosch (2002; based on Monte 
Carlo realizations of EPS), respectively. 
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Fig. 11. — Predicted and simulated median halo MAHs for sim- 
ulations LCDMl-3. Symbols and lines arc the same as in Figure 
10. 

for fow-mass halos in some CDM models, they have very 
different asymptotic behaviors at high redshift, because 
the W02 model is an exponential function of z while ours 
is a power-law function of z. 

Several months after this current paper was submitted 
to the Journal and the Internet, McBride et al. (2009) 
modeled the distribution of halo MAHs of the Millen- 
nium simulation by fitting individual MAHs with a two- 
parameter function, instead of the one-parameter func- 
tion of W02. They found that our model prediction for 
the mass accretion rate to have a slightly steeper z depen- 
dence than their Equation (9) where our median value is 
within 20% of their median mass accretion rate at z ~ 
but exceeds theirs by a factor of ~ 2 at z ~ 4. The 
reason for this large discrepancy is unknown, however, 
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Fig. 12. — Predicted and simulated median halo MAHs for sim- 
ulations SCDMl-2 and OCDMl. Symbols and lines are the same 
as in Figure 10. 




z+1 

Fig. 13.— Predicted and the PINOCCHIO median halo MAHs 
for a cosmology model the same as LCDMl— 3. Black symbols 
represent the median redshifts of the halo MAHs output by the 
PINOCCHIO code automatically, while solid red lines are our 
model predictions. Four different PINOCCHIO runs with box sizes 
of 25, 100, 300, and 1200?i-iMpc and particle numbers of 256^, 
256^, 512^, and 256"^ are analyzed here. The vertical dotted line 
indicates the redshift 4. 

it is worthy to point out that since we combine a set 
of high-resolution simulations, the dynamical ranges ex- 
plored here are larger than the Millennium Simulation 
and our halo samples are well enough for analyzing the 
median value. To test our model further, we utilize the 
Lagrangian semianalytic code PINOCCHIO,^ proposed 
by Monaco et al. (2002) for identifying dark matter ha- 
los in a given numerical realization of the initial den- 
sity field in a hierarchical universe. Figure 13 compares 
our model predictions with the median redshifts of the 
MAHs produced by the PINOCCHIO code, for a cos- 

® A public version of PINOCCHIO is available from the Web 
site: http:/ /adlibitum. oats. inaf.it/monaco/Homepage/Pinocchio. 
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mology model the same as LCDMl-3 and very similar to 
that of the Millennium Simulation. Here four different 
PINOCCHIO runs with box sizes of 25, 100, 300, and 
1200/i-iMpc and particle numbers of 256^, 256^, 512^, 
and 256^ are analyzed and, to avoid artifact, only those 
MAHs which converge among runs of different boxsizes 
are plotted. Again, at all redshifts available, our model 
predictions match the median values of the PINOCCHIO 
MAHs very well in halo mass range probed by the Millen- 
nium Simulation. The fluctuations on the median MAHs 
in the smallest box are due to the small number statis- 
tics. Note that the PINOCCHIO MAHs arc automati- 
cally output by the PINOCCHIO code without any post- 
processing and so are much independent of our numerical 
simulations and of our definitions of halo merger trees. 

It is interesting to examine the scaling relations ob- 
tained above more closely. At sufficiently high z when 
lg(5c(z) - lg(5c(zobs) > 0.272/w(zobs), Equation (10) 
gives p = 0, and Equation (12) reduces to a simpler 
form. In this case, w(z,M) = (5c(z)/s(Af) — 5.85 cor- 
responds to dlgcr(M)/dlg(5c(2) = 1, or equivalently to 
d\%v I d\%8\[z) — 0. The scaled mass, y^, is then in- 
dependent of time, so is w in the SF case. This cor- 
responds to the high-redshift asymptotic behavior dis- 
cussed in the last subsection. For halos more massive 
than this characteristic scale, the model predicts that 
d\'ga(M')l dX^bciz) > 1 and that ^/i/ increases with time. 
For scale- free cases, this means that the most massive ha- 
los will grow faster and faster, as demonstrated in the two 
upper panels of Figure 10, while for ACDM and OCDM 
cases, this kind of rapidly growing halos are very rare be- 
cause neither the 5c~l + z relation nor the a-M relation 
is a simple power law. Both the decrease of cosmologi- 
cal density parameter Q„i with time and the increase of 
power index of the power spectrum with mass scale will 
slow down the halo mass growth rate, as shown above. 
This is why at the present time halos on nearly all mass 
scales show two distinct phases in their mass growth his- 
tories, as found by Z03b and Z03a. 

For the SF models, we can obtain some useful analytic 
formulae, because here a{M) has a simple power-law de- 
pendence on M: a oc M~*^"+'^)/^. First, using the high- 
z asymptote, w{z,M) = 5.85, we have 6c{z)/a{M) = 
5.85 X l0dis'y/d\g^ ^ 5 35 ^ io-("+3)/6 and 1^ = 5.85^ x 
j^Q-("+3)/3_ 'pjjis means that the mass of a median pro- 
genitor at sufficiently high redshift is a fixed fraction of 
the characteristic mass: M = 5.85^/("+3)(M*/10). For 
n — 1, 0, —1 and —2, these progenitors have a typical 
mass M/M, = 1.415,3.422,20.02, and 4008.07, respec- 
tively. Second, because d\ga/d\gm = —{n + 3)/6 is a 
constant in an SF model, we have dig s{M) / d\gdc{z) = 
d\ga{M)/dlgSc{z) = w{z,M)/5M = dc{z)/s{M)/5M 
when p{z) = 0. The solution is s{M) = ddz) /5.85 + C 
with C a constant along each median history. This so- 
lution describes the MAHs in the SF models quite well, 
as shown in Figure 8. The s{M)-Sc{z) relations in the 
CDM models are steeper than that given by this simple 
model, as shown in Figure 9, because the effective power 
spectrum index, which comes into the definition of s{M) 
increases with the mass scale. 

Before ending this section, we would like to point out 
that the scaling relations obtained here may also apply 
for individual halos. As shown in Figure 14, the linear re- 




FlG. 14. — Simulated mass accretion histories for individual halos 
in simulations LCDMl— 3. For clarity, histories of different final 
masses are presented with different colors. 



lation between s{M) and Sc{z) is also a good approxima- 
tion for individual halos, although the slope may change 
from halo to halo. Thus, one may model the individ- 
ual MAHs using a straight line with its slope changing 
from halo to halo. The scatter in the MAHs for halos of 
a given mass may then be modeled through the distri- 
bution of the slope. Since the linear relation is a good 
approximation for different halo masses and in different 
cosmologies, as shown in Figures 8 and 9, this kind of 
modeling is expected to be valid for different halos and 
in different cosmological models. We will come back to 
this point in a forthcoming paper. 

4. A UNIVERSAL MODEL FOR THE DARK MATTER 
HALO CONCENTRATIONS 

4.1. The Mass and Redshift Dependence of Halo 
Concentration 

Median values of the halo concentrations measured for 
the CDM models are plotted as symbols in Figures 15- 
17. The results obtained here from simulations LCDMl- 
3 are in excellent agreement with the results obtained 
in Z03b. In these plots, we also compare our simula- 
tion results with the three much-used empirical mod- 
els. The first one is the NEW model, which relates the 
halo characteristic density at Rg to the universe den- 
sity at the time when 50% of the halo mass is already 
in progenitors of 1% of the halo mass or bigger. In 
agreement with BOl, our results show that the NEW 
model not only fails to predict correctly the redshift de- 
pendence of halo concentration, but also fails to pre- 
dict the concentration at z = accurately. The model 
of EOl matches our results better, especially at z = 1 
and z = 2, but it also fails to match the concentration- 
mass relation, particularly at high redshift. The model of 
BOl has several versions of model parameters (see Web 
site http://www.physics.uci.edu/~bullock/CVIR/). We 
adopt the original version with F = 0.01 and K = 3.75, 
where F and K are parameters in Equations (9) and (12) 



Dark Halo Growth and Structure 



11 



z = 0.00 z=1.03 




— Bullock 2001 (K=3.75) \ \ 

Eke 3001 ^- 
This work 



1 ll[ l l ll| I [ l l ll l l l I lllllll| I lllllll| I lllllll| I I lll l lll l I I l llllll| I lllllll| I lllllll| l ' 
fi _ z = 2.0B z = 3.94 




10"° 10" 10'2 10'= 10'« 10" 10" 10'= 10" lO'* lO'S 
M (h-' Mg) M (h-' MJ 

Fig. 15. — Predicted and simulated median concentration as a 
function of halo mass for simulations LCDMl-3. Black symbols 
are results measured in simulations and error bars the standard 
deviation among halos of the same mass divided by the square 
root of the halo number in the bin, and these data are in excellent 
agreement with those presented in Figure 2 of Zhao et al. (2003b). 
Solid red lines are our model predictions, while short-dashed green, 
long-dashed dark blue, and dot-dashed blue lines are predictions of 
the NFW model, the BOl model, and the EOl model, respectively. 
The BOl model has several versions of model parameters and here 
we adopt the version with F = 0.01 and K = 3.75. 

of their paper. This version matches well our simula- 
tion results for LCDMl-3 at 2: = 0. However, it fails to 
match the concentration-mass relation for massive halos, 
especially at high z. If F = 0.01 and K = 3.4 is used, 
as suggested for total halo population lately by Maccio 
et al. (2007) and James Bullock on his Web site, the 
model underpredicts the concentration by an amount of 
10% even for low-mass halos at z = 0. The conflicts be- 
tween early model predictions and simulation data have 
already been discussed by Z03b, and were also confirmed 
recently by Gao et al. (2008) who used the Millennium 
Simulation to carry out an analysis very similar to that 
in Z03b. As mentioned in the introduction, Gao et al. 
tried some revisions of the NFW model and found that 
the revised version still fails to match their simulation re- 
sults. The simulation results presented here are in good 
agreement with the results in Z03b and in Gao et al. 
(2008). Since we combine a set of high- resolution simu- 
lations, the dynamical ranges explored in Z03b and here 
are larger than the Millennium Simulation. With our 
new simulation data of the SCDM model, we find that 
the BOl cannot even match our data at z = unless pa- 
rameters in their model are adjusted. As in the LCDM 
case, the BOl model also fails to account for the redshift 
dependence of halo concentration in the SCDM model. 

Z03a found that the scale radius Rg is tightly corre- 
lated with the mass Mg enclosed by it. With this tight 
correlation, one can predict halo concentration from the 
MAHs. Using the MAH model given in last section, we 
can predict the halo concentration according to the Z03a 

^ Although in the published paper of BOl K = Ais adopted, their 
latest Web site has deleted all description about this value and 
claim that K = 3.75 should be the original version corresponding 
to the published paper. 
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Fig. 16. — Predicted and simulated median concentration as a 
function of halo mass for simulations LCDM4-6. Symbols and error 
bars are the same as in Figure 15. Solid red lines are our model 
predictions, while short-dashed green, long-dashed dark blue, and 
dot-dashed blue lines are predictions of the model of Maccio et al. 
(2007), the BOl model, and the Z03a model, respectively. 
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Fig. 17. — Predicted and simulated median concentration as a 
function of halo mass for simulations SCDMl-3. Symbols and lines 
are the same as in Figure 15. 

model. The prediction of this model is shown in Figure 
16 as the dot-dashed line. For comparison, we also show 
predictions of the BOl model and the model of Maccio 
et al. (2007) in the same figure. As one can see, the 
prediction of the Z03a model is more accurate than the 
other two models, particularly for high mass halos. 

Although the model of Z03a matches well the simula- 
tion data, it is not easy to implement. Here we use our 
simulation results to find a new model that is more ac- 
curate and easier to implement. Let us first consider an 
SF case with a given linear power index. In this case, the 
halos of a given mass have a uniquely determined time 
tp when their main progenitors reach a fixed fraction 
e.g., 50%, of their final mass. There is thus a one-to-one 
correspondence between the MAH and tp- If the con- 
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Fig. 18. — Concentration of the final off'spring halo as a function 
of the universe age when its main progenitor first gained 4% of 
the final halo mass. Results are shown for scale-free cases with 
n = 0, —1, —2 and for LCDM3. Different colors represent different 



snapshots, and the red line is for z ■■ 
by Equation (13). 



: 0. The solid curves are given 



centration of a halo is fully determined by its MAH, one 
may build a model to link tp and the final halo concen- 
tration. In principle, such a relation can be found for tp 
at an arbitrarily fixed fraction of the final halo mass, but 
the resulting relation may be very different for scale-free 
models of different power indices. What we are seeking 
is a value of F with which the relations between c and 
tp are the same for all the scale-free models. 

After many trials, we find that the concentration of 
a halo is tightly correlated with the time io.o4 when its 
main progenitor gained 4% of its final mass. This relation 
is almost identical for n = 0, —1 and —2, as shown in 
Figure 18, but is slightly lower for the n = 1 model. 
Since for realistic cosmological models, the effective slope 
of the linear power spectrum on scales of halos that can 
form is always less than 0, we exclude the result of the 
n = 1 model in modeling the c-to.o4 relation. We find 
that the relation given by the other three models can be 
accurately described by the following simple expression: 



[4« + (t/toMf 



4 X [1 + (t/3.75io.o4) 



8.411/8 

Nontrivially, this same relation also applies very well to 
the CDM models, as shown in the small window of Figure 
18. 

From our MAH model described in Section 3, we can 
easily compute the time to.o4 for a halo of mass M. We 
can then predict its concentration straightforwardly by 
using Equation (13). The concentration-halo mass rela- 
tions so obtained are shown as the solid lines in Figure 
19 for the SF models, and in Figures 15-17 for the CDM 
models. Comparing the model predictions with the cor- 
responding simulation data, we see that our model works 
accurately for all the models and for halos at different 
redshifts. The typical error of our prediction is less than 
5%. Somewhat larger deviations (~ 10%) seen for halos 
of mass ^ 10^^ Mq at redshift z ~ 3 could be a result of 
numerical inaccuracy in the simulations, because these 
halos are not well relaxed. Compared with early mod- 
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Fig. 19. — Halo median concentration as a function of halo mass 
for scale-free models with n = 1, 0, —1, and —2. Colored symbols 
are results measured in simulations and error bars are the standard 
deviation among halos of the same mass divided by the square 
root of the halo number in the bin. Smooth lines are our model 
predictions. For clarity, both the simulation data and the model 
predictions have been shifted in horizontal direction by -1, 0, 1.5, 
and 3.5 for SFl, SF2, SF3, and SF4, respectively. 

els, ours is clearly a very significant improvement. In 
Appendix B, we describe step-by-step how to use Equa- 
tion 13 to predict the concentration-mass relation at any 
given redshift for a given cosmological model. 

It should be pointed out that both our model and the 
model of NEW are based on the correlation between halo 
concentration and a characteristic formation time. How- 
ever, the two models have several important differences. 
First, in the NEW model the formation time was defined 
as the epoch when half of the halo mass M has been 
assembled in its progenitors of masses exceeding O.OIM, 
while in our model the time is defined as the epoch when 
its main progenitor has gained 4% of the halo mass. Sec- 
ond, the ways to relate the concentration to the charac- 
teristic time are also different. NEW assumed that the 
inner density at Rg of a halo is proportional to the mean 
density of the universe at the formation time, while we re- 
late the halo concentration and the characteristic time by 
Equation (13). Finally, NEW used the extended Press- 
Schechter formula to compute the formation time, while 
we use our model of MAHs to calculate the time to.o4- 
These differences make a very big difference in the model 
predictions, as shown above. 

4.2. The Evolution of Halo Structure Along the Main 

Branch 

The model described above can also be used to predict 
how halo structural properties, such as c, Rs, Ps, Ms, 
and Vs, and hence density profile evolve along the main 
branch. For a given MAH, M{z), we can estimate to.04 
for the current halo, M{z) [corresponding to the virial 
radius R{z) and the circular velocity 1^(2;)], at the red- 
shift z in question, and use the model described above to 
obtain c at z. Since the virial radius R is determined by 
M{z), we can then obtain Rs, Ps, Mg, and Vg through 
Rs = R/c, Equations (4) and (5), respectively. The solid 
lines in Figure 20 show the model predictions in compar- 
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Fig. 20. — Predicted (smooth lines) and simulated (symbols con- 
nected with lines) evolution of the median c, ps, Ms, and Vs 
for halos of nine different final masses in simulations LCDMl-3. 
The increasing direction of final halo mass among different lines in 
left/right panels is downward/ upward. Here simulated ps, Ms, and 
Vs are estimated with the simulated median c and the predicted 
halo mass. For comparison, evolution of circular velocity at the 
virial radius, V , is also plotted in the lower right panel as dashed 
lines. 

ison with the simulation results of LCDMl-3 shown by 
circles connected by lines. Clearly, our model also works 
very well in describing these evolutions. These results 
demonstrate that, for a given halo, our model cannot 
only predict how its total mass grows with time, but can 
also predict how its inner structure changes with time. 
Thus, one can plot its NFW density profile at any point 
of the main branch, and can obtain the density evolution 
in spherical shell of any radius. 

As one can see, for low-mass halos, ps, Ms, Vg, hence 
also Rs, all remain more or less constant at low redshift, 
suggesting that the inner structures of these halos change 
only little in the late stages of their evolution. Since by 
definition the virial radius R increases with time, the 
concentration of such halos increases rapidly with de- 
creasing redshift. Note that the circular velocity at the 
virial radius actually decreases with time for such halos 
at late time, as shown by the dashed lines in the lower- 
right panel of Figure 20. On the other hand, for massive 
halos, Ps, Ms, Vs, and Rs all change rapidly with red- 
shift even at z ^ 0, implying that the inner structures of 
these halos are still being adjusted by the mass accretion. 
Therefore, the concentration of such halos increases very 
slowly or even remains constant with redshift, much dif- 
ferent from the W02 model which argues that c oc a. All 
these behaviors had been illustrated clearly in Z03a and 
Z03b. As discussed in Z03a, these different behaviors are 
mainly due to the fact that massive halos are still in their 
early growth phases, which is characterized by rapid halo 
growth dominated by major mergers, effectively reconfig- 
uring and deepening the gravitational potential wells and 
causing the collisionless dark matter particles to undergo 
dynamical relaxation and to mix up sufficiently to form 
the core structure, while low-mass halos have reached the 
late growth phase, which is characterized by slower quies- 
cent growth predominantly through accretion of material 
onto the halo outskirt, little affecting the inner structure 
and potential. 

Figure 21 presents the model predictions of the evolu- 
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Fig. 21. — Predicted (smooth lines) and simulated (symbols con- 
nected with lines) evolution of the median c for halos of different 
final masses in simulations SFl-4. The increasing direction of final 
halo mass among different lines is from top to bottom. For com- 
parison, assumption of W02, c oc 1/(1 + z), is also plotted in each 
panel as the dashed line. 

tion of median c in comparison with the simulation re- 
sults of SFl-4. It is very interesting that the simulated 
median concentrations of massive halos in each SF model 
also remain constant but the values are not 4 except for 
SF4 (n = —2). This is not surprising: as shown in Figure 
10, the mass accretion of SF4 halos, which is the fastest 
in these SF models, is very similar to the CDM cases and 
thus triggers dynamical relaxation and core-structure 
formation, resulting a constant concentration of about 
4; while for the rest SF models, mass accretion is more 
or less slower and even the massive straight-line-MAH 
halos is in their slow growth regime (t/^o.o4 > 3.75; see 
Figure 18) and, according to Equation (13), progenitors 
on these straight line MAHs, which have time-invariant 
V^o.o4j will all have constant but n-dependant concentra- 
tions. According to Section 3.1, all median MAHs in SF 
cases asymptote M (x (1 -I- 2:)"^/^^+"^ and, in an Einstein 
de Sitter universe, universe age t oc (1 + z)^^'^. Com- 
bining these two relations with Equation (13), we find 
halo median concentrations also have a unique asymp- 
totic behavior at high redshift: c ^ 29.37, 12.61, 5.477, 
and 4.007 for n = 1, 0, —1, and —2, respectively. These 
numbers are well verified in Figure 21 (except that in the 
SFl case c should be shifted a little bit when compared 
with the simulation results because a short-wavelength 
cutoff in the linear power spectrum have been included in 
this simulation) and, for comparison, assumption of W02, 
coc l/(l-|-z), is also plotted in each panel of the figure. 
Actually, for this kind of self-similar straight-line MAHs, 
any concentration model which claim that halo concen- 
tration is tightly correlated to its MAH should predict 
a constant concentration, as required by logic. For the 
most massive and so very rare halos which grow faster 
and faster, as demonstrated in the two upper panels of 
Figure 10, our model predict a decrease of concentration 
with time, as its i/to.o4 is diminishing. This very in- 
teresting behavior is again supported by the simulation 
results and worth further detailed study. 

4.3. Comparison with the Zhao et al. (2003a) Model 
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As mentioned above, the Z03a model for the hafo con- 
centrations is based on the tight correlation between Rs 
and Ms- In Figure 22 we show the Ms-Rs relation ob- 
tained from the simulations used in this paper. Clearly, 
for a given mass, Mg and Rs are tightly correlated, and 
the relation is well described by Ms oc i?^" with a ~ 0.55. 
The value of a is between the values 0.48 and 0.64 ob- 
tained in Z03a for the late slow and early fast growth 
regimes, respectively. This suggests that the model we 
arc proposing here is closely related to that of Z03a. To 
show this more clearly, let us start with Equation (10) in 
Z03a: 



[ln(l + c) 


-c/(U 


- c)]c-3« 
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[ln(l + Co) - 
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(14) 

The function of c in square brackets on the left side be- 
haves like a power law for moderate concentrations and 
then we will use a power law function instead. Suppose 
that we have a halo whose main branch reaches 4% of 
its current mass at a time to.o4 and assume that the con- 
centration of the progenitor at t = toM is 4. Substitut- 
ing quantities of subscript in the above equation with 
quantities at i = torn gives 



P0.04 
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obtain 



0.603 for moderate c < 25. We can then 
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assuming a — 0.55. This is very similar to relation (13) 
at t > 3.75io.o4i indicating that the model proposed here 
is consistent with that of Z03a. Indeed, for a given MAH, 
the redshift when t — 3.75 x io.o4 is uniquely deter- 
mined. We found that this redshift separates well the 
fast growth regime {t < 3.75to.o4) from the slow growth 
regime {t > 3.75io.o4)- In the fast growth regime all halos 
have about the same median concentration c ~ 4 inde- 
pendent of halo mass, while in the slow growth regime 
the concentration scales with time as ^ 4(t/3.75io.o4)^ °^- 
This is exactly the proposal made in Z03a, Z03b that con- 
centration evolution of a halo depends much on its mass 
accretion rate and the faster the mass grows, the slower 
the concentration increases. 

5. CONCLUSIONS AND DISCUSSION 

A large amount of observations have constrained cos- 
mological parameters and the initial density fluctuation 
spectrum to a very high accuracy. However, cosmological 
parameters change with time and the power index of the 
power spectrum varies with mass scale dramatically in 
the so-called concordance ACDM cosmology. Thus, any 
successful model for its structural evolution should work 
well simultaneously for various cosmological models and 
different power spectra. 

In this paper, we use a large set of A^-body simula- 
tions of various structure formation models to study the 
MAHs, the mass and redshift dependence of concentra- 
tions and the concentration evolution histories of dark 
matter halos. We find that our simulation results cannot 
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Fig. 22. — Inner mass Ms vs. inner radius Rs along the simulated 
median MAH obtained from LCDMl-3. Results are shown for 
halos with nine different final masses, the same as those shown in 
Figure 20. 



be described by the much-used empirical models in the 
literature. Using our simulation results, we developed 
new empirical models for both the MAHs and the con- 
centration evolution histories of dark matter halos, and 
the latter can also be used to predict the mass and red- 
shift dependence of halo concentrations. These models 
are universal, in the sense that the same set of model pa- 
rameters works well for different cosmological models and 
for halos of different masses at different redshifts. Our 
model predictions also match the PINOCCHIO MAHs 
very well, which are automatically output by the La- 
grangian semianalytic code PINOCCHIO without any 
postprocessing and so are much independent of our nu- 
merical simulations and of our definitions of halo merger 
trees. These models are also simple and easy to imple- 
ment, making them very useful in modeling the growth 
and structure of the dark matter halo population. 

We found that, in describing the median MAH of dark 
matter halos, some degree of universality can be achieved 
by using Sc{z), the critical overdensity for spherical col- 
lapse, to represent the time, and cr(M), the rms of the 
linear density field on mass scale M, to represent the 
mass. This is consistent with the Press-Schechter (PS) 
theory, in which the dependence on cosmology and power 
spectrum is all through Sc{z) and a{M). A more univer- 
sal relation can be found by taking into account the local 
slope of the power spectrum on the mass scale in consid- 
eration. This is again expected, because the growth of an 
average halo in a given cosmological background is de- 
termined by the perturbation power spectrum on larger 
scales. Indeed, in the extended PS theory, the average 
mass accretion rate of a halo of a given mass depends 
both on (j(M) and its derivative. However, the exact 
dependence on the shape of the power spectrum in our 
empirical model is actually different from that given by 
the extended PS theory, as demonstrated clearly by the 
discrepancy between our results and the model of van den 
Bosch (2002), which is based on the extended PS theory 
(see Figures 10-12). We find that another modification 
is required in order to correct the deviation of halos cor- 
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responding to relatively low peaks. We suggest that such 
deviation is due to the neglect of the large-scale tidal field 
on the formation of dark matter halos, as found in Wang 
et al. (2007). Thus, the empirical model we find for the 
halo accretion histories may help us to develop new mod- 
els for the mass function and merger trees of dark matter 
halos. In addition, we find that a linear relation between 
s{M) [defined in Equation (8)] and Sdz) is also a good 
approximation for individual halos, although the slope 
may change from halo to halo. Thus, the MAHs of indi- 
vidual halos may be modeled with a set of straight lines, 
in the s{M)-Sc{z) space, with different slopes specified 
by a distribution function. Research along these lines is 
clearly worthwhile, and we plan to come back to these 
problems in forthcoming papers. 

Our other finding is that the concentration of a halo 
is strongly correlated with to.o4, the age of the universe 
when its progenitor on the main branch first reaches 4% 
of its current mass. This is consistent with the general 
idea that the structure of a dark matter halo is correlated 
with its formation history. The concentration of a halo 
selected at a given redshift is determined by the ratio of 
the universe age at this redshift, t, to to.o4- If this ratio 
is smaller than 3.75, the halo is in the fast growth regime 
and its concentration is expected to be c ~ 4, indepen- 
dent of its mass. For a halo in the fast growth regime, 
its first 4% mass settled down at a late epoch (of age 
> t/3.75) and so has low density. In addition, its current 
accretion rate is large in the sense that it only takes a 
time interval of much less than 2.75t/3.75 to gain the 96 
percent of mass remaining. Such accretion is expected to 
have a significant impact on the inner structure that is 
not very dense, and so the inner structure of this halo is 
adjusted constantly by the mass accretion. In contrast, 
for a halo in the slow growth regime (f/to.04 > 3.75), its 
first 4% mass settled down at an early epoch when the 
universe age was smaller than f/3.75, and so this fraction 
of the mass is expected to settle into a very dense clump. 
Also, for such a halo, the current accretion rate is small 
in the sense that it takes a time interval of more than 
2.75f/3.75 to gain the 96% of mass remaining. Such slow 
accretion is not expected to have a significant impact on 
the inner dense structure that has already formed. All 
these are in qualitative good agreement with the results 
presented in ZG3a. Our model can also be used to predict 
accurately how halo structural properties, such as c, Rg, 
Ps, Ms and Vs, and hence its density profile evolve along 
the main branch, making it very useful for our under- 
standing of the formation and evolution of dark matter 
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halos. Assuming that the NFW profile is a good approxi- 
mation for the halo density profile, the model can predict 
too the MAHs, the mass and redshift dependence of con- 
centrations and the concentration evolution histories for 
different halo definitions. In a forthcoming paper, we will 
apply this model to individual halos. 

We found that in a self-similar cosmology halos on 
some characteristic mass scale have a straight-line me- 
dian MAH and a time-invariant median concentration. 
All halos asymptote these median features at high red- 
shift: halos under this mass scale will grow slower and 
slower on average and their concentration will increase 
with time; while halos above this mass scale will grow 
faster and faster on average and their concentration will 
decrease. On the other hand, in ACDM and OCDM 
cases, both the decrease of cosmological density parame- 
ter ^Irn with time and the increase of power index of the 
power spectrum with mass scale will slow down the halo 
mass growth rate, inducing that at the present time ha- 
los on nearly all mass scales show two distinct phases in 
their mass growth histories, as found by Z03b and Z03a. 
Z03a and Z03b also pointed out that when halo mass ac- 
cretion rate reaches some criteria, the fast mass growth 
will trigger dynamical relaxation and core-structure for- 
mation, resulting in a constant concentration of about 
4. Here we found another mechanism that will lead to a 
time-invariant concentration along the main branch of a 
halo merger tree: a self-similar straight- line MAH will in- 
duce a constant concentration, as required by logic, and 
when the mass accretion rate is below the above criteria, 
the value of this constant is determined by the accretion 
rate and is always larger than 4. 
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APPENDIX 

In the appendices A and B, we give a detailed user's guide on how to use the results of the present paper to compute 
the median mass accretion history and the concentration of dark matter halos for a given cosmological model. 

CALCULATING THE MEDIAN MASS ACCRETION HISTORY OF HALOS 

For a given cosmological model, one can calculate Sc{z) = 6[.[^lm{z),^A{z)]/ D{z) for a given redshift z, where 
S[.[Q,miz),flA{z)] is the conventional spherical collapse threshold at redshift z, D{z) is the linear growth factor nor- 
malized to 1 at redshift z — Q, and 5c{z) is the collapse threshold for the density field linearly evolved to ^; = 0. Given 
a linear power spectrum at 2; = 0, one can also calculate the rms density fluctuation a{m) for a spherical volume of 
mass m, as well as ,s(m). Once these quantities arc obtained, one can calculate the median mass growth history of 
dark matter halos of mass Mobs at redshift Zobsj i-S-, the median mass M(2:|2:obs) -^obs) of their main progenitors at a 
higher redshift z > Zohs- The step-by-step procedure is as follows. 

1. With a{Mohs)-, s(Mobs) and 5c(-2obs) for halos of mass Mobs at redshift Zo\ys, one can calculate w{zohs-, Mohs) 
according to Equation (7) and the shift p(-2obs, -^obs, -^obs) according to Equation (11). Let z = Zohs and M = Mobs 
temporarily. 

2. One can obtain the accretion rate dlgcr(M)/(ilg(5c(-z) at this redshift by substituting the above w and p into 
Equation (12). 

3. At an incrementally higher redshift z' = z + Az, one can calculate the now critical collapse threshold 5c{z') at 
z' , and its change A\g5c{z) between z' and z, and hence obtain the change in \ga{M) through Alg(7(M) = 
d\ga{M)/d\g5c{z) x A\g5c{z). One can then get a new (t(M'), with M' being the median mass of the main 
progenitors at z' . The mass M' can be obtained from a through the function a{m.). 

4. In order to trace the history further backward, one needs to calculate w{z', M') according to Equation (7) and 
calculate the shift p{z' , Zobs, Mobs) according to Equation (10). Let z = z', M = M', and repeat steps 2 and 3. 

Then step by step, one can trace the MAH backwards to high redshift. One can also calculate the halo radius R{z) 
and halo circular velocity V{z) along this MAH according to Equations (2 and 3) by setting Ah = Avir as defined in 
Section 1. 

CALCULATING THE HALO CONCENTRATION 

For a given MAH ending at [zobs^Mobs] in a given cosmology, one can predict concentrations for all the main 
progenitors at different redshifts in addition to the final offspring, i.e., one can obtain the concentrations along the 
MAH, c(2;|zobs, -^obs)- Then many other useful quantities can be calculated as shown below. There are two steps in 
computing the concentrations. 

1. Prom a given MAH, M(z|zobs) -^obs)) one can compute Zo,q4,, the redshift at which the main progenitor of the 
main progenitor halo at z{> Zohs) grows to 0.04M, where M is the mass of the main progenitor halo at z. 

2. In the given cosmology, one can calculate the ages of the universe, t and to. 04; at these two redshifts, z and Zo.o4; 
and obtain c(z|zobs, Mobs) by substituting them into Equation (13). 

3. Furthermore, one can calculate the characteristic inner quantities of the halo along the MAH, such as 
i?s(z|zobs, Mobs), Ps{z\zohs,Mohs)-, Ms(z|zobs, Mobs) and T4(z|zobs, Mobs), according to the definition of c and 
Equations (2-5), and hence can plot the NFW density profile at any point of the MAH according to Equation 

(!)• 

In the literature, there are several definitions for halo radius. Assuming that the NFW profile is a good ap- 
proximation for the halo density profile, one can predict the MAHs and the concentration evolution histories for 
different halo definitions, because halos of different definitions have the same inner structures (siich as Rs, Ps and 
Ms) although their boundaries are different. For example, according to Equations (4 and 5), one can calculate 
C200c(-2|-2obs,-M^obs), Af200c('Z|-2obs, Afobs) , V2ooc{z\zobs, Mobs) for halo definition Ah{z) = 2QQ/Vtm{z) and one can also 
calculate C2oom(2;|2;obs, Afobs), M2oorn{z\zohs-,Mohs), V200m(-2|^obs, M)bs) for halo definition A/,(z) = 200. 
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Note that c{z\zobs, Afobs) obtained with steps 1 and 2 is the concentration of a main progenitor of a halo of mass Afobs 
at Zobs- As we pointed out in Section 3.2, this population of main progenitors of mass M at z is not statistically the 
same as the whole population of halos of mass M at z, although in some cases the difference between them is negligible. 
The former population resides in environment of more frequent mergers and so has slightly smaller concentration. In 
any case, one can easily compute c(M, z) for the whole population of halos of mass M at ^; if one sets ^obs = z and 
Mobs = M in the above calculation. 

A CALCULATOR, A USER-FRIENDLY CODE AND TABLES AVAILABLE ON INTERNET 

A calculator which allows one to interactively generate the median MAHs, the concentration evolution 

histories and the mass and redshift dependence of concentrations of dark matter halos is provided at 
http : //www. shao . ac. cn/dhzhao/mandc.html, together with a user-friendly code to make the relevant calculations. 
The calculator and the code can be used for almost all current cosmological models with or without inclusion of 
the effects of baryons in the initial power spectrum. We also provide tables for median concentrations of halos with 
different masses at different redshifts in several popular cosmological models. Detailed instructions on how to use the 
calculator, the code and the tables can be found at the Web site. 



